/****************************************************************************
 * nbody.c
 * An implementation of the naive O(n^2)/timestep numerical integration
 * approach to the N-body problem. 
 ****************************************************************************/

#include "nbody.h"

/****************************************************************************
 * Allocate and setup initial conditions for the bodies
 ****************************************************************************/
Body *init_bodies(unsigned int num_bodies, int init_type) {
  int i;
  double n = num_bodies;
  Body *bodies = (Body *) malloc(num_bodies * sizeof(Body));
  
  for (i = 0; i < num_bodies; i++) {
    switch (init_type) {
    case INIT_LINEAR:
      bodies[i].mass = 1.0;
      bodies[i].pos[X] = i / n;
      bodies[i].pos[Y] = i / n;
      bodies[i].vel[X] = 0.0;
      bodies[i].vel[Y] = 0.0;
      break;
    case INIT_SPIRAL:
      bodies[i].mass = (n - i) / n;
      bodies[i].pos[X] = (1 + i / n) * cos(2 * M_PI * i / n) / 2;
      bodies[i].pos[Y] = (1 + i / n) * sin(2 * M_PI * i / n) / 2;
      bodies[i].vel[X] = 0.0;
      bodies[i].vel[Y] = 0.0;
      break;
    }
  }

  return bodies;
}

/****************************************************************************
 * Verify that the simulation is running correctly:
 * it should satisfy the invariant of conservation of momentum
 ****************************************************************************/
int check_simulation(Body *bodies, int num_bodies) {
  int i, check_ok;
  double momentum[2] = {0.0, 0.0};

  for (i=0; i<num_bodies; i++) {
    momentum[X] += bodies[i].mass * bodies[i].vel[X];
    momentum[Y] += bodies[i].mass * bodies[i].vel[Y];
  }

  check_ok = ((abs(momentum[X]) < EPSILON) && (abs(momentum[Y]) < EPSILON));
  if (!check_ok) printf("Warning: total momentum = (%3.3f, %3.3f)\n", momentum[X], momentum[Y]);
  return check_ok;
}

/****************************************************************************/

int main(int argc, char **argv) {
  struct timeval start_time, end_time;
  double time;

  int num_bodies, num_steps;
  int i, j, k;
  double dt = 1.0, dv[2];
  double r[2], dist, force_len, force_ij[2], tot_force_i[2];
  Body *bodies;

  /* Parse command line */
  const char Usage[] = "Usage: nbody <num bodies> <num_steps>\n";
  if (argc < 2) {
    fprintf(stderr, Usage);
    exit(1);
  }

  num_bodies = atoi(argv[1]);
  num_steps = atoi(argv[2]);

  /**************************************************************************
   * Initialize bodies and allocate internal data structures
   *************************************************************************/

  //printf("Initializing with %d bodies and %d time steps\n", num_bodies, num_steps);

  bodies = init_bodies(num_bodies, INIT_SPIRAL);
  check_simulation(bodies, num_bodies);

  /* SEQUENTIAL EXECUTION */

  /* Start timing */
  //printf("Running "); fflush(stdout);
  gettimeofday (&start_time, 0);

  /**************************************************************************
   * Main loop
   *************************************************************************/

/* Store in force_ij[] the force on body i by body j */
#define Calc_Force_ij() \
  r[X] = bodies[j].pos[X] - bodies[i].pos[X]; \
  r[Y] = bodies[j].pos[Y] - bodies[i].pos[Y]; \
  dist = r[X] * r[X] + r[Y] * r[Y]; \
  force_len = GRAV_CONST * bodies[i].mass * bodies[j].mass / (dist * sqrt(dist)); \
  force_ij[X] = force_len * r[X]; \
  force_ij[Y] = force_len * r[Y]

/* Update velocity and position of body i, by numerical integration */
#define Step_Body_i() \
  dv[X] = dt * tot_force_i[X] / bodies[i].mass; \
  dv[Y] = dt * tot_force_i[Y] / bodies[i].mass; \
  bodies[i].pos[X] += dt * ( bodies[i].vel[X] + dv[X] / 2 ); \
  bodies[i].pos[Y] += dt * ( bodies[i].vel[Y] + dv[Y] / 2 ); \
  bodies[i].vel[X] += dt * dv[X]; \
  bodies[i].vel[Y] += dt * dv[Y]

  for (k = 0; k < num_steps; k++) {
    //printf(". "); fflush(stdout);

    {

    for (i = 0;  i <num_bodies; i++) {
      tot_force_i[X] = 0.0;
      tot_force_i[Y] = 0.0;
      /* Compute total force f(i) on each body i */
      for (j = 0; j < num_bodies; j++) {
	    if (j == i) continue;
	    Calc_Force_ij();
	    tot_force_i[X] += force_ij[X];
	    tot_force_i[Y] += force_ij[Y];
      }
      Step_Body_i();
    }

    }
    
  }

  /* Stop timing */
  gettimeofday (&end_time, 0);
  time = end_time.tv_sec - start_time.tv_sec;
  time = time + (end_time.tv_usec - start_time.tv_usec) / 1e6;

  /* END SEQUENTIAL EXECUTION */
  
  check_simulation(bodies, num_bodies);

  //printf("done!  interaction rate: \n%6.0f\n", num_bodies * (num_bodies - 1) * num_steps / time / 1000);
  
  printf ("%d:%f\n", 0, time);
  
  return(0);
}

